The interplay between shell effects and electron correlations in quantum dots 
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We use the Path Integral Monte Carlo method to investi- 
gate the interplay between shell effects and electron correla- 
tions in single quantum dots with up to 12 electrons. By use 
of an energy estimator based on the hypervirial theorem of 
Hirschfelder we study the energy contributions of different in- 
teraction terms in detail. We discuss under which conditions 
the total spin of the electrons is given by Hund's rule, and the 
temperature dependence of the crystallization effects. 

PACS numbers: 73.20.Dx, 71.45.Lr, 75.30.Fv 73.23-b 



I. INTRODUCTION 

The advances in nanofabrication of the last years 
opened the goal to build 2D quantum dots (QDs) and 
quantum dot molecules (QDMs) - artificial mesoscopic 
semiconductor structures of selectable shape and size - 
as containers for a controllable fixed number of elec- 
trons Recently, depending on the strength and 
shape of the effective confining potential, the formation 
of spin density waves (SDWs) an d Wigner crystals 
[pLpl in QDs and QDMs has been predicted by differ- 
ent groups with different theoretical approaches. Hirosc 
and Wingreen || argue that SDWs are reproducible arte- 
facts of spin density functional calculations. For a 2D 
parabolic confining potential the accordance of the spin- 
configuration with Hund's Rule has been predicted by 
Koskinen, Manninen, and Reimann Q and questioned by 
Yannouleas and Landman g]. All these effects are gov- 
erned by the intriguing interplay between shell effects, 
the pure coulomb repulsion, and the fcrmionic repulsion 
due to the Pauli exclusion principle and depend strongly 
on the values of the interaction parameters in the com- 
monly for single QDs assumed Hamiltonian 
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where k is the dielectric constant, m* is the effective 
mass, and loq defines the strength of the confining po- 
tential. 

Apart from the interesting physical questions that arise 
for quantum dots the reliable prediction of their proper- 
ties is an ultimate test of modern methods in quantum 
chemistry. Due to the compared to atoms very shallow 
confining potential long range electron interactions and 
correlations play an important role in QDs and QDMs. 
Therefore it is misleading to name them artificial atoms 



and molecules. Well established and very elaborate meth- 
ods of quantum chemistry might fail in describing them 
properly. Hartree-Fock and spin density functional meth- 
ods use single Slater determinants or sums of them to ap- 
proximate the many-body wave function. In spin-density 
functional methods the approximation of the functional 
for the exchange correlation energy adds another 

source of uncertainty and systematic errors to this ap- 
proach. The Path Integral Monte Carlo method (PIMC) 
used in this paper samples the full many-body wave func- 
tion instead. 

In contrast to density functional methods with PIMC it 
is possible to study the temperature dependent properties 
of QDs. The reason why PIMC is not yet a standard 
method of quantum chemistry is its numerical limitation 
due to the fermion sign problem. The rapidly increasing 
power of modern computers resizes this limitation. In 
sec. H we briefly summarize our implementation of PIMC 
and comment on how to limit the numerical deficiencies 
due to the fermion sign problem. 

We apply PIMC to calculate the electron density and 
two-particle correlation functions for quantum dots with 
up to 12 electrons. To compare to various experimen- 
tal as well as to other theoretical studies we use differ- 
ent dielectric constants k and strengths of the confining 
parabolic potentials. The calculated addition energies 
are in very good agreement with the experimental find- 
ings of Tarucha et al. ||. 

For N — 6 we investigate the temperature dependence 
of the Wigner crystallization (WC). 



II. NUMERICAL METHOD 

For a system of N electrons with position eigenket 
I %i,Si) (sj = ±i for spin-up and spin-down electrons) 
in an external potential the Feynman path integral can 
be written as [9 11 
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and the boundary condition Xj(M + 1) = Xj(l). M is 
the number of so-called timeslices of the Feynman paths. 
In the limit M — > oo Eq. (||) becomes exact. For quan- 
tum dots the space dimension is d—2 and the (2NM)- 
dimensional integral given in (|^) can be evaluated by 
standard Metropolis Monte Carlo techniques. Due to the 
determinant the integrand is not always positive and the 
expectation value of an observable X(x.) depending only 
on position operators has to be calculated using 
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where X g is the value of the observable X and W g is 



the value of the integrand in (|2j) in the g'th Monte Carlo 
step. Eq. (|]) reveals a severe problem connected with 
the path integral for fermions which is commonly de- 
noted as the fermion sign problem (see e.g. p^-[T^] ) . It 
can be shown that the ratio between integrands with pos- 
itive sign (VF + ) and negative sign (W~) is approximately 
given by @Jl5| 
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where Ep and Eb are the ground state energies of the 
Fermi system and the corresponding Bose system. It is 
now obvious that the statistical error in (||) grows rapidly 
for small temperatures T . Moreover the energy difference 
(Ep — Eb) will grow with increasing system size causing 
an increase of the statistical error. 

Within PIMC the calculation of the kinetic energy ex- 
pectation value is another critical task. This is merely 
due to the fact that the Monte Carlo calculation is usu- 
ally done in position space and that the discretization of 
the paths allows a number of different approaches to cal- 
culate the expectation value of a momentum dependent 
operator. A number of various different energy estima- 
tors has been discussed in the past 

To avoid these difficulties we developed a procedure 
which allows the calculation of all energy expectation val- 
ues from the knowledge of the pair correlation functions 



= (S(r- | x.-xj I)) 
and the radial density functions per electron 
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where q is the probability of finding electron i in distance 
r from the center. 

Due to the particle symmetry we have 
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Utilizing the hypervirial theorem of Hirschfelder JTgj ] the 
energy can be written as a sum of ten parts p3] 



E = E L + E L + E ll + E tL + E iL 

+ E lot + E iot + E lot + E pot + E pot 



(10) 
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While in density functional approaches the calculation of 
the kinetic energy and the exchange correlation energy 
is a major topic and subject to permanent discussion, 
within the path integral approach these energies are in- 
cluded in a natural way. 

However, the systematic error arising from the limited 
number of timeslices M and the statistical error of the 
Monte Carlo calculation have to be controlled carefully. 
We checked our algorithm extensively using eight non- 
interacting fermions in a parabolic trap as a test system. 
We found that at low temperatures where the ratio of 
signs is around 0.99, convergence can only be achieved 
obeying the following rules: 1) The determinants have to 
be calculated very accurately using a more costly algo- 
rithm with pivoting. 2) The completely uncorrelated gen- 
eration of the Monte Carlo steps is essential, i.e. the coor- 
dinate to be moved should be chosen randomly. Moving 
the particle coordinates using always the same sequence 
produces inaccurate results. 3) A good random number 
generator with a completely uncorrelated sequence in all 
significant bits of a 64 bit real number should be ap- 
plied. We therefore developed a 53 bit random number 
of Marsaglia-Zaman type instead of using one of the 
standard 24 or 32 bit random number generators coming 
with standard system libraries. 4) Further, to improve 
the convergence a number of different Monte Carlo steps 
can be applied, i.e. moving single time slices, moving 
complete particle paths and parts of a path. 
Our Fortran-code is completely parallelized using MPI 
and Lapack. 
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III. RESULTS 

To compare our PIMC calculations to experimental 
data we calculated the addition energies 

AE = -E/v+i — 2£jv + -EW— l (H) 

of a QD with up to 11 electrons using the material con- 
stants m* — 0.067m and k = 12.9 for GaAs as given by 
Hirose and Wingreen || . It is assumed that these param- 
eters mimic the experimental setup of Tarucha et al. ||] 
reasonably well. The strength of the harmonic potential 
is fixed at ?uvo = 3.0 meV. The resulting effective atomic 
units are E* H = 10.955 meV for the Hartree energy and 
<2q = 10.1886 nm for the Bohr radius. The Boltzmann 
constant is k B = 7.8661 x 10~ 3 £^/K. 

We performed PIMC simulations for quantum dots 
with different spin configurations at a fixed temperature 
of 10 K. Due to the fermion sign problem the number of 
Monte Carlo steps necessary to push the statistical error 
of the total energy, which has been calculated properly 
from 25 uncorrelated subsequences of MC steps, into the 
range of 0.1 percent is extremely high. The number of 
Monte Carlo steps ranged between 2.5 billion steps per 
particle coordinate for N < 6 and about 10 billion steps 
for N = 12. Fig. [I] displays the addition energies for 
quantum dots with up to 11 electrons. The circles indi- 
cate the results from our path integral calculations at 10 
K, the squares are results of spin density functional calcu- 
lations of Hirose and Wingreen || , and the triangles are 
the experimental results of Tarucha et al. Both theo- 
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FIG. 1. Addition energies for quantum dots with up to 
11 electrons. The circles indicate the results from our path 
integral calculations at 10 K, the squares are results of spin 
density functional calculations of Hirose and Wingreen, and 
the triangles are the experimental results of Tarucha et al. 
The error bars for PIMC would be of the size of the solid 
circles and are therefore omitted. 

retical calculations reproduce the general ^-dependence 
of the addition energies in great detail. Tarucha et al. 



give an estimate of the electron temperature in their ex- 
periments of T = 0.2 K. For computational reasons our 
PIMC calculations are performed at 10 K and it cannot 
be expected that the absolute energy values agree as well 
as the K DFT calculations with the experimental re- 
sults. However, it should be noted that PIMC correctly 
predicts the drop in the addition energy from N = 7 to 
N ~ 8 while the DFT calculations fail at this point. 

The inset in Fig. [I] displays the total spins of the spin 
configurations with lowest energy as found in DFT and 
PIMC at 10 K. In DFT calculations (0 K) the spin config- 
uration of the ground state is determined by Hund's rule 
for up to 22 electrons. In contrast, in our PIMC calcu- 
lation at 10 K the total spin is not always in accordance 
with Hund's rule. For TV = 4 we checked the tempera- 
ture dependence of the spin configuration. At 5 K the 
energy of the spin configuration is 0.01 E\ higher than 
the spin 1 energy indicating a temperature dependence 
of the favored spin configuration. 

As an important fact we note that the iV-dependence 
of the addition energies is not affected by the actual spin 
configuration. The situation is quite similar to that in 
transition metal clusters with extreme small energy dif- 
ferences between states with significantly different mag- 
netic moments [ p2[ . 

As can be inferred from Fig. |(a) the radial spin den- 
sities are significantly different for both spin configura- 
tions. The total potential energy for the spin 1 configu- 
ration is about 0.07 meV lower than that of the spin 
configuration. At 10 K this is overcompensated by an 
0.27 meV higher kinetic energy (see Tab. ||). Although 
the kinetic and potential energies for different total spins 
significantly differ the total energies are almost equal. 
Similar situations are found for larger N. 

For convenience and easy comparision we determined 
the value of the dimensionless density parameter r s , 
which is sometimes used to characterize quantum dots 
(see e.g. f§) to be r 5 =4.19 for N = 4. 

The energies given in Tab. | correspond to the integrals 
in Eq. [l0| The total kinetic energy which is the sum of 
all E£ in terms is always positive while some of the ad- 
dends might be negative. Tab. Q reveals that larger total 
spins result in larger kinetic energies. The total potential 
energy is almost unchanged, the larger contribution from 
the trap potential is compensated by a smaller contribu- 
tion from the Coulomb repulsion. We note that the ratio 
£w between the kinetic energy and the total energy is 
considerably larger for N = 4 than for N — 9 reflecting 
the looser binding of the smaller system. 

Next we consider the dependence of the Wigner crys- 
tallization on the temperature and the choice of the mate- 
rial constants. The localization of the electrons in space is 
commonly referred as Wigner crystallization. For quan- 
tum dots the occurrence of well separated humps in the 
radial electron density and the pair correlation functions 
has been interpreted as WC. However, it is a nontrivial 
task to find a general parameter identifying if an electron 
system is crystallized or not. From a solid state physics 
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iV T =2,N l =2 


N T =3,N l =l 


N J =5,N l =4 


N T =6,N l =3 


E to t 


40.83 


41.03 


169.25 


169.82 


Skin 


7.33 


7.60 


18.97 


19.57 


-Epot 


33.50 


33.43 


150.28 


150.26 


Ew 


0.18 


0.19 


0.11 


0.12 


E, ■ 

km 


8.03 (8.03) 


12.66 (8.44) 


35.44 (14.18) 


43.74 (14.58) 


1-1 1 

Et 

kin 


8.03 (8.03) 


3.55 (7.10) 


27.30 (13.65) 


19.39 (12.93) 


Ei . 

km 


-1.31 (-2.62) 


-3.97 (-2.65) 


-11.20 (-2.24) 


-16.87 (-2.25) 


^kin 


-6.11 (-3.05) 


-4.64 (-3.09) 


-25.84 (-2.58) 


-23.31 (-2.59) 


Et- 

km 


-1.31 (-2.62) 


0.00 (0.00) 


-6.72 (-2.24) 


-3.39 (-2.26) 


Epot 


33.50 


33.43 


150.28 


150.26 


Epot 


8.03 (4.01) 


12.66 (4.22) 


35.44 (7.09) 


43.74 (7.30) 


Epot 


8.03 (4.01) 


3.55 (3.55) 


27.30 (6.82) 


19.39 (6.46) 


E n 

^pot 


2.62 (2.62) 


7.94 (2.65) 


22.40 (2.24) 


33.73 (2.25) 


^pot 


12.21 (3.05) 


9.28 (3.09) 


51.69 (2.58) 


46.62 (2.59) 


Epot 


2.62 (2.62) 


0.00 (0.00) 


13.44 (2.24) 


6.78 (2.26) 



TABLE I. Kinetic and potential energies as well as 
Eyf=Ekin/ Et t in meV for different electron configurations 
at 10 K (the numbers in brackets are the single particle ener- 
gies). 



eters. Fig. |3j displays the radial densities and pair corre- 
lation functions for six electrons with S=0, hu=5 meV 
and k=3.0, 6.0, and 12.9. Of course, for stronger elec- 
tron repulsions (small k) the electron distributions are 
broadened. The qualitative picture of the distributions 
is merely the same. For all k shell effects indicated by 
off-center maximums of the radial density occur. How- 
ever, only for k—3 and 6 we observe a maximum at r=0. 
From our point of view it cannot definitely be decided 
from this figure if a system is Wigner crystallized or not. 
As a parameter reflecting the strength of the WC we em- 
ploy the ratios between the kinetic and the total energies 
Ew = E kin /E tot which are 0.07, 0.10 and 0.14 for «=3.0, 
6.0 and 12.9. Although the radial distribution function 
for k =12.9 is quite narrow, the relative mobility for the 
electrons indicated by Ew is twice as large as for k =3. 
The underlying physical process can be understood intu- 
itively. Due to the stronger electron-electron repulsions 
the electrons are fixed in an energetical favorable geomet- 
ric configuration and as a consequence thereof the relative 
kinetic energy is reduced and the difference between the 
pair correlation functions of equal and oppposite spin al- 
most vanishes (see Fig. . It is an interesting and to our 
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FIG. 2. Radial density per electron (a) and pair correlation 
functions (b) for 4 electrons and total spins S—0 and 5=1. 
The material constants are k = 12.9 and Tiuj — 3.0 meV. 



point of view the electrons should have a low mobility, i.e. 
a small kinetic energy, and should not interchange their 
lattice positions. For fermions the localization of single 
electrons does not make any sense, and, as stated above, 
even the decomposition of the many body wave function 
in sums of determinants of single particle wave functions 
is probably a too rough approximation for QDs. These 
facts limit the analogies between crystallization in solids 
and electron systems and make the term crystallization 
itself somehow misleading. We therefore view Wigner 
crystals as states of the many body wave function with 
a relatively low kinetic energy. 

First we consider the strength of the Wigner crystal- 
lization depending on the choice of the interaction param- 




b) 



5 10 15 20 25 30 

r[aj 



10 20 30 40 50 60 

r [a *] 



FIG. 3. Radial density per electron (left) and pair corre- 
lation functions (right) for N = 6, S = 0, hcoo=5 meV and 
dielectric constants k=3.Q, 6.0, and 12.9 at 10 K. The radial 
density for k=12.9 is scaled by a factor of 3 and the pair 
correlation functions are scaled to have a maximum value of 
±1. 

knowledge open question, if the crystallization of elec- 
trons can be viewed as a phase transition. We therefore 
consider next the temperature dependent properties of a 
quantum dot with N=6, S=0, k=3 and hu—5 meV. The 
results for temperatures between 10 and 150 K are pre- 
sented in Tab. n[ Most notably Ew increases relatively 
smoothly from 0.07 at 10 K to 0.22 at 150 K. Within 
our numerical accuracy the caloric curve does not show 
any evidence of a phase transition. The transition from a 
crystallized state to an electron fluid seems to be squashy. 
Of course, from our calculations we cannot exclude that 
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FIG. 4. Radial density (left) and pair correlation functions 
(right) for N = 6, S = 0, hu!o=5 meV, k — 3.0, and temper- 
atures T = 10, 30, 60, and 150 K. 



A comparision to other QMC methods seems to be in 
order here. Even our most complicated simulations took 
less than two hours on a Cray T3E with 62 processors. 
Taking the advantages of modern computer power and 
optimized software, the brute force PIMC applied here 
is able to produce very precise results. Different algo- 
rithms have been published to improve the performance 
of path integral methods. A recent one is the Multi-Level 
Blocking method published by Mak et al. Q . Simula- 
tions using this algorithm are expected to converge better 
than our direct treatment since the fermion sign problem 
is avoided partially. The published results of the treat- 
ment of quantum dots by Egger et al. || do not confirm 
this expectation. Obviously the advantages of the new 
method are more than compensated by other numerical 
problems, maybe due to the energy estimator used p4[ . 
Nevertheless, a combination of the Multi-Level-Blocking 
method and our technique might result in a very powerful 
tool. 



a phase transition exists for larger TV or different inter- 
action parameters. Fig. ^ displays the radial electron 
densities and the total electron pair correlation function 
for different temperatures. Up to 60 K the radial density 
shows clear geometric structure effects with two maxi- 
mums while at 150 K only a smooth curve resembling to 
a simple gaussian remains. 
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N T =3, N l =3 



10 K 
30 K 
60 K 
90 K 



T = 120 K 
T — 150 K 



Ew 
0.07 
0.08 
0.12 
0.16 
0.19 
0.22 



g kin [meV] 
17.70 
22.93 
35.39 
49.44 
64.10 
78.49 



ffpot [meV] 
246.53 
249.53 
257.67 
267.01 
276.84 
286.88 



fftot [meV] 
264.32 
272.46 
293.06 
316.46 
340.93 
365.35 



TABLE II. Kinetic and potential energies for different tem- 
peratures and spin configuration iV — 3. The Hartree 
energy is E^ — 202.558 meV, k — 3 and hu>o = 5.0 meV. Ew 
is the ratio between the kinetic and the total energy. 



IV. CONCLUSION 

In conclusion, we have found that despite of the noto- 
rious fermion sign problem PIMC is capable to answer 
interesting questions for strongly correlated electron sys- 
tems like QDs and QDMs. For QDs PIMC reproduces 
correctly the experimental addition energies. Our tem- 
perature dependent calculations give new insights into 
the process of WC. For the 2-dimcnsional QDs a ratio 
Ew = Eki n /E toi below 0.1 seems to indicate WC both 
for k and temperature dependent calculations. However, 
reagarding this aspect a more firm classification parame- 
ter e.g. similary to the Lindemann criterion is desirable. 



[1] N. Zhitenev, R. Ashoori, L. Pfeiffer, and K. West, Phys. 

Rev. Lett. 79, 2308 (1997). 
[2] S. Tarucha, D. Austing, and T. Honda, Phys. Rev. Lett. 

77, 3613 (1996). 
[3] C. Yannouleas and U. Landman, Phys. Rev. Lett. 82, 

5325 (1999). 

[4] M. Koskinen, M. Manninen, and S. Reimann, Phys. Rev. 

Lett. 79, 1389 (1997). 
[5] R. Egger, W. Hausler, C. Mak, and H. Grabert, Phys. 

Rev. Lett. 82, 3320 (1999). 
[6] K. Hirose and N. Wingreen, Phys. Rev. B 59, 4604 

(1999). 

[7] D. Pfannkuche, V. Gudmundsson, and P. Maksym, Phys. 

Rev. B 47, 2244 (1993). 
[8] R. Pino, Phys. Rev. B 58, 4644 (1998). 
[9] M. Takahashi and M. Imada, J. Phys. Soc. Jpn. 53, 936 

(1984). 

[10] M. Takahashi and M. Imada, J. Phys. Soc. Jpn. 53, 3765 
(1984). 

[11] P. Borrmann and E. Hilf, Z. Phys. D 26, S350 (1993). 

[12] D. Ceperley, Phys. Rev. Lett. 69, 331 (1992). 

[13] A. Lyubartsev and P. Vorontsov-Velayaminov, Phys. 

Rev. A 48, 4075 (1993). 
[14] I. Morgenstern, Z. Phys. B 77, 267 (1989). 
[15] W. H. Newman and A. Kuki, J. Chem. Phys. 96, 356 

(1992). 



[16] J. Cao and B. Berne, J. Chem. Phys. 91, 6359 (1989). 
[17] A. Giansanti and G. Jacucci, J. Chem. Phys. 89, 7454 
(1988). 

[18] H. Kono, A. Takasaka, and S. Lin, J. Chem. Phys. 88, 
6390 (1988). 

[19] J. O. Hirschfelder, J. Chem. Phys. 33, 1462 (1960). 
[20] H. Heinze, P. Borrmann, H. Stamerjohanns, and E. Hilf, 

Z. Phys. D 40, 190 (1997). 
[21] P. Borrmann, in preparation (unpublished). 
[22] K. Lee and J. Callaway, Phys. Rev. B 48, 15358 (1993) 
[23] C. Mak, R. Egger, and H. Weber-Gottschick, Phys. Rev. 

Lett. 81, 4533 (1998). 
[24] D. Ceperley, Phys. Rev. Lett. 69, 331 (1992). 



6 



